Signatures of shear thinning— thickening transition in dense athermal shear flows 
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In non-equilibrium molecular dynamics simulations of dense athermal shear flows, we observe the 
transition from shear thinning to shear thickening at a crossover shear rate 7 C . Shear thickening 
occurs when d ," n . a , > 2 with T g the granular temperature. At the transition, the pair distribution 
function shows the strongest anisotropy. Meanwhile, the dynamics undergo apparent changes, sig- 
nified by distinct scaling behaviors of the mean squared displacement and relaxation time on both 
sides of 7 C . These features serve as robust signatures of the shear thinning— thickening transition. 

PACS numbers: 83.60.Rs,83.60.Fg,83.10.Rs,83.80.Fg 



Complex fluids such as colloids, gels, emulsions, foams, 
and granular materials exhibit intriguing and compli- 
cated rheological phenomena subject to shear. Unlike 
newtonian fluids whose shear viscosity is independent of 
the shear rate, complex fluids can behave shear thinning 
or thickening under proper conditions, characterized by 
the decrease or increase of the viscosity with increas- 
ing the shear rate [l|-|4|- Shear thinning and thickening 
are attractive research topics with important applications 
[2, la] • In practice, we always wish paints to shear thin, 
while shear thickening is desired in the design and man- 
ufacture of smart materials such as soft body armors. 

The underlying mechanisms of shear thinning and 
thickening have been debated over decades. Early simu- 
lations have suggested that the layering of particles along 
the direction of the shear flow contributes to shear thin- 
ning [fj, lZ[, which has been suspected as an artifact of 
profile biased thermostat |8J |9( and been challenged by 
most recent studies (lfj, [ll|. A recent measure of the 
microscopic single-particle dynamics has proposed that 
shear thinning results from the decrease of the entropic 
force contribution [3|. Besides, shear thinning can be 
boosted by the presence of yield stress which on the con- 
trary impedes the emergence of shear thickening [12J . 
Compared to shear thinning, our understanding of shear 
thickening is even poorer, because fewer systems exhibit 
shear thickening and shear thickening usually happens 
at high shear rates where probing is difficult. Multi- 
ple mechanisms, e.g. the order-disorder transition and 
formation of hydroclusters [2|, |7|, ll3J], have been pro- 
posed to explain shear thickening. Recent studies have 
also unveiled possible links between shear thickening and 
jamming of constituent particles due to dilation [14l4l7|. 



Shear thickening has been observed in both Brownian 
and Non-Brownian suspensions [3j, |l2|, [l7( and in simu- 
lations with and without taking hydrodynamics into ac- 
count [2, lD, lfj, LL8(- It then remains an open question 
whether shear thickening originates from the same mech- 
anism for various systems. 



In this letter, we study the rheology of planar shear 
flows of athermal granular systems via non-equilibrium 
molecular dynamics simulations. By applying the sim- 
ple model without the interference of hydrodynamics and 
thermostat, we observe the transition from shear thin- 
ning to shear thickening. This transition is accompanied 
with some robust signatures which may not be limited to 
athermal systems and may thus provide us with a general 
picture of shear thickening. 

Our systems are two-dimensional L x L squares con- 
sisting of disks with an identical mass m. A half of the 
disks have a diameter of <r, while the other half have 
(Tl = 1.4cr. Lees-Edwards boundary conditions [19( are 
applied with the shear being imposed in the a;— direction. 
SLLOD equations of motion assuming a linear velocity 
profile [20( are employed: 
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where fi = (xi,y,i) and Vi = {v X i,v y i) are the location 
and random velocity of particle i, 7 = -3 is the shear 
rate with 7 the shear strain, F(j is the force acting on 
particle i by particle j, and the sum is over all parti- 
cles j interacting with particle i. The force F%j includes 



two parts, the elastic force F? 



-VVij and damping 



force F^j = -£m{Vij + A fy. lJ x) [2l|, where Vi 



< 13 , ^j , 



and 

i/ij are the interaction potential, relative random veloc- 
ity, and y— distance between particles i and j, and £ is the 
damping coefficient. Wceks-Chandler-Anderson (WCA) 
potential, i.e. repulsive Lennard- Jones potential, is ap- 



plied: V i3 
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when the sepa- 



ration between particles i and j, r^, is smaller than the 
sum of their radii 07,, and zero otherwise. We use a, 
m, and e as the units of length, mass, and energy. The 
units of temperature and time are e/ks and a/y/e/m, 



respectively, where fcs is the Boltzmann constant. 

We apply Gear predictor-corrector algorithm to in- 
tegrate Eqs. ([J} and © at a constant packing frac- 
tion (j) = 0.85 above the T — jamming transition at 
4> c f=a 0.84 [22M25J I . The systems at rest are jammed solids 
obtained from L-BFGS energy minimization |26j. Data 
are collected and averaged over time after the systems 
have been sheared over a time much longer than the char- 
acteristic time scale 1/7 and steady shear flows without 
the memory of their initial states are achieved. The shear 
stress T, xy is calculated from 
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where Y? xy and T, x x y 



are the ideal gas stress and excess 
stress from particle interactions. Correspondingly, the 
shear viscosity can be divided into two parts: 



f] = f] ld + rf x 



K d y/i + ^x X y/i- 



(4) 



We vary the number of particles TV from 256 to 4096 to 
verify that our results do not show significant system size 
dependence. Here we only show results for TV = 1024. 

In a steady shear flow, energy injected into the flow by 
the shear force is dissipated by the damping force, from 
which we obtain 
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FIG. 1: (color online) (a) Viscosity 77 calculated from Eq. (|4} 
against the shear rate 7 at f = 0.1 (black circles), 0.2 (red 
squares), 0.5 (blue diamonds), and 1.0 (green triangles). The 
solid curves are eddy viscosity rj e calculated from Eq. ([6]). 
(b) Shear rate dependence of the ideal gas viscosity rf d (blue 
squares) and excess viscosity r) ex (red circles) at £ = 0.1. The 
black curve is the total viscosity, (c) Shear rate dependence of 
the four components of the eddy viscosity, rj T ^ (red circles), 
— r l_ 2 (blue squares), 1 97 x | (green diamonds), and r/ (ma- 
genta triangles), at £ = 0.1. The solid curve is the total eddy 
viscosity, (d) Shear rate dependence of the rate of granular 
temperature increase A at the same damping coefficients as 
panel (a). Inset: Comparison of r\ and A for a thermostated 
system at -J2_ £. »& = 0.1. 



where the sum is over all interacting particle pairs. Sim- 
ple calculations of Eq. ([5]) lead to the eddy viscosity [8J 
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On the right hand side of Eq. ([7]), 2% is the average co- 



Si=i ^f denotes the 



ordination number, and T g = 2J J3I 
granular temperature, i.e. effective temperature at short 
time scales [2l|, [27j . When introducing T g in Eq. (0 , 
we assume that the random part of the kinetic energy, 



,iV 



yii—i ^mv 2 , is equally partitioned to all particles. 



The decomposition of the viscosity by Eqs. (jj) and 
© enables us to find out the actual source of the shear 
stress in control of the rheology. Equation (|4j is the mi- 
croscopic expression of the viscosity from virial theorem, 
while Eq. ([6]) is straightforwardly derived from energy 
balance. In homogeneous steady flows we would expect 
that r\ = r/ e . 

Figure[lja) shows flow curves at four different damping 
coefficients ranging from 0.1 to 1. There is a crossover 
shear rate, j c , separating shear thinning at low shear 
rates from shear thickening at high shear rates. Our sys- 
tems at rest have a nonzero yield stress. Shear thinning 
is consequently expected at low shear rates [12j. When 
scaling 7 by £, we notice that % ~ £. Scaling collapse 
at high shear rates in the shear thickening regime is also 
observed when rj/£ is plotted against 7/^, which approx- 
imately obeys a power law scaling 77 ~ j 2 . This scaling 
implies that the shear thickening flows here are not Bag- 
noldian (77 ~ 7) [171, |28[. To check if our results are due 
to the Lees-Edwards boundary conditions, we have also 
studied systems confined between two walls at constant 
volume with the top wall moving in the x— direction at a 



constant speed and observed similar results. As will be 
discussed in the following, our shear thickening flows are 
associated with high temperature and high pressure gas- 
like states induced by very high shear rates, which are 
in the same regime as previous simulations [8|, [9( but in 
different regime from recently reported Bagnoldian shear 
thickening flows arising from unjamming— jamming tran- 
sition at low shear rates [17j . 

As compared in Fig. 03a), rj « r\ e as expected, ex- 
cept in the vicinity of j c where Arj = r/ e — rj > 0. A?y 
tends to increase with increasing £. We have verified that 
A 77 > is not a transient behavior by observing no time 
evolution. In the shear rate regime where Arj is appar- 
ently greater than 0, the kinetics show visible anisotropy, 
e -9- J2i v xi > Y^i v yi- A s will be discussed, there exists 
structure anisotropy in the same regime as well. The 
anisotropy should thus account for the nonzero A77. Al- 
though r\ and r\ e are not completely equal, they exhibit 
the same j c , so this inequality does not affect our dis- 
cussions about the shear thinning— thickening transition 
using Eqs. (j4j and (J6J respectively. 

Let us first see what Eq. Q tells us. In Fig. QJb), we 
compare rf d with rf x . At low shear rates, the random 
motion of particles is so slow that tf x ^> rf ■ With in- 
creasing the shear rate, rf d grows up and eventually beats 
rf x . jf d = rf x at a crossover shear rate 7 9 . When 7 > j g , 
the system behaves effectively as a high temperature and 
high pressure gas with £ _1 and 7" 1 the dominant time 
scales, which may be the cause of the high shear rate 
scaling collapse shown in Fig. [Ha). Because 7 g > j c , 
Eq. (j4]) does not give us any clue about how the shear 
thinning— thickening transition happens. 

In contrast, Eq. (JU) is significantly useful to reveal the 
possible source of shear thickening. In Fig. QJc) , we show 
all the four viscosity components in Eq. (jB]). In the whole 
shear rate regime studied here, rj Q and rj_ 1 are small and 
negligible, r/ T 2 > 0, and r\_ 2 < 0. At low shear rates, 
\tj_ 2 \ is smaller than but comparable to tj t 2 , and shows 
similar shear rate dependence to rf^ 2 . Near j c , however, 
rf[ 2 becomes significantly larger than \rj_ \. Shear thick- 
ening is thus mainly determined by rj r 2 . Assuming that 
the average coordination number zj, does not vary largely 
with the shear rate, which is actually true for our sys- 
tems, Eq. (O implies that if the granular temperature T g 
varies faster than 7 2 shear thickening would occur. 

In athermal shear flows with a constant damping co- 
efficient, T g increases with increasing the shear rate. To 
estimate the rate of T g increase, we define a quantity 



A 



d(lnT a ) 
d(ln7) 



and plot it against the shear rate in Fig. [TJd) . 
At low shear rates, A w 1.5, indicating that T g ~ 7 15 . 
This scaling law stops working near j c . When 7 > j c , 
A > 2 and shear thickening happens, exactly as expected 
from our analysis of Eq. (JB]). 

We have now established a direct link between the 
granular temperature and shear thickening for our model 
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FIG. 2: (color online) (a) — (c) Snapshots at £ = 0.1 and 
7 = 0.01, 0.1, and 1 with -y c « 0.1. (d)-(f) Pair distri- 
bution functions g(x,y) for systems shown in (a) — (c). The 
value of g(x,y) is quantified by the color, (g) — (h) Pair dis- 
tribution functions parallel and perpendicular to y = x, g" (r) 
and g ± (r), at 7 = 0.01 (black solid), 0.1 (red dashed), and 1 
(blue dot-dashed), (i) Anisotropy of g(f) measured by a at 
£ = 0.1 (black circles), 0.2 (red squares), 0.5 (blue diamonds), 
and 1.0 (green triangles). 



systems. Is this link general or only specific to our sys- 
tems? Does it depend on particle interactions? First 
of all, we have verified that systems with harmonic, 
Hertzian, and Lcnnard- Jones interactions come to the 
same conclusion. Next, let us see if A > 2 is also the con- 
dition of shear thickening in other typical model systems. 
Besides our athermal model, typical models to study 
shear flows include the foam model [29|, |30( , Langevin dy- 
3l|,[32|, and thermostated shear flows (a. \8Lffl, |27| . 



namics 

For the foam model described by the equation of motion, 
(vi = J2j F-tj j where ( is equivalent to £m in our model, 
simple calculations result in rj e ~ X^i ^li~ 2 ■ We can find 
at once that shear thickening happens when dffn ' )"' > 
2, equivalent to A > 2 except that particle inertia is ig- 
nored. When inertia is not neglected, A > 2 can be 
straightforwardly obtained as well. For Langevin dy- 
namics, a random force Ri(t) is added to the equation of 
motion for the foam model. Because Vi and Ri are uncor- 
rected |32|, we would expect the same result. For ther- 
mostated shear flows, we have already known that when 
T g is fixed shear thickening cannot happen [6[ , no matter 
if it is an artifact of thermostat [19( . Even with more re- 
alistic thermostat, it has been reported that T g is much 
higher than the bath temperature [19(. These previous 
observations may hint that shear thickening is related to 



the granular temperature in thermostatcd systems. In 
the inset to Fig. QJd), we show both 77 and A against 
the shear rate for a thermostated shear flow governed by 
!u — m Sj F ij ~~ i v yi% ~ av yiVi where a maintains a 
constant kinetic energy in the y— direction [27f . For this 
model, energy balance results in r\ e ~ a 12i v yi A f^ 2 ^ fr° m 
which A > 2 cannot be recognized at all. Interestingly, 
A > 2 still signifies shear thickening. 

For foam model and Langevin dynamics, A > 2 is a 
straightforward consequence of energy balance for shear 
thickening to occur. For our athermal model, this conse- 
quence is not so obvious because the precondition is that 
•q T is significantly larger than the other three terms in 
Eq. ©, which turns out to be true. For thermostated 
shear flows, A > 2 is totally unpredictable. We thus pro- 
pose that A > 2 is a generic signature of shear thickening 
at high shear rates. 

The granular temperature is difficult to measure ex- 
perimentally. It is then interesting to know if there are 
any experimentally accessible signatures associated with 
the shear thinning— thickening transition. As illustrated 
in panels (a) — (c) of Fig. [2] the system undergoes appar- 
ent changes from shear thinning to shear thickening. In 
the shear thinning regime, the system is roughly uniform 
in space, while in the shear thickening regime particles 
tend to form instantaneous clusters with large overlaps 
and leave behind more voids. We are thus inspired to 
search for possible structural or dynamical signatures of 
the shear thinning— thickening transition. 

Recent studies have attempted to build the 
brid ge between shear rheology and microstructure 
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L 1 



33 35]. 



Here we measure the pair distribu- 
function for large particles, g(f) = g(x, y) = 

XW/zFy2iY,jj:iS('r(x,y)-fi j )y with special atten- 
tion to its anisotropy, where the sums are over all 
large particles and (.) denotes the time average. The 
contour plots of g(r) shown in panels (d) — (f) of Fig. [2] 
indicate that g{f) is asymmetric in the vicinity of the 
shear thinning— thickening transition: the contours are 
elongated in the direction of y = x and compressed in 
the direction perpendicular to y = x. 

In order to quantify the asymmetry, we show in pan- 
els (g) and (h) of Fig. [2] the pair distribution function 
parallel and perpendicular to y = x, denoted as g"(r) 
and g ± (r), respectively. The first peak of g ± (r) is higher 



than that of y"(r). 

.11 



Moreover, 



..II 



C r's , where r~" and 
r's denotes the length at which g ± (r) and g"(r) start to 
be greater than zero. These features imply that when 
the structure is asymmetric particles tend to have larger 
and more uniform overlaps in the direction perpendic- 
ular to y = x. We define a = rl/rjr to characterize 
the asymmetry. Figure [5Ji) shows that with the increase 
of shear rate a first increases and then drops after a 
maximum. The maximum emerges approximately at j c . 
Shear thinning— thickening transition is thus signified by 
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FIG. 3: (color online) (a) Mean squared displacement (Ay 2 ) 
in the shear thinning (magenta solid symbols) and thickening 
(green empty symbols) regimes at £ = 0.1. The black circles 
are at 7 ~ j c . The shear rate increases from the right to the 
left in the inset. In the main panel, t is scaled by the relax- 
ation time r. (b) Shear rate dependence of the relaxation time 
r at £ = 0.1 (black circles), 0.2 (red squares), 0.5 (blue dia- 
monds), and 1.0 (green triangles). The solid lines show power 
law behaviors in the shear thinning and thickening regimes. 



the most pronounced structure anisotropy. 

In the inset to Fig. [3ta), we show the y— component of 
the mean squared displacement (MSD) for large particles, 
(Ay 2 ), where (.) denotes the particle and ensemble aver- 
age. There is a short time ballistic motion ((Ay 2 ) ~ i 2 ) 
followed by diffusion ((Ay 2 ) ~ t) at long times. Ap- 
parently, the relaxation time r decreases with increasing 
the shear rate. When we plot the MSD's against t/r, 
where r is determined from (Ay 2 (r)) rj 10er 2 , as shown 
in Fig. |3|a), the MSD's in the shear thinning regime col- 
lapse nicely onto the same curve, while the ballistic parts 
of the shear thickening curves are still apart. The scaling 
collapse of shear thinning curves implies that particles 
move ballistically to approximately the same distance. 
In the shear thickening regime, however, particles move 
ballistically to a longer distance with increasing the shear 
rate. As shown in Fig.^c), particles form instantaneous 
clusters in shear thickening flows and leave lots of voids. 
If clustered particles move collectively at short times, the 
voids allow the particles to move ballistically to a longer 
distance, which may be the picture of the short time dy- 
namics of shear thickening flows. 

The distinct dynamics between shear thinning and 
thickening flows can also be identified from the shear rate 
evolution of the relaxation time. As shown in Fig. [2Jb), 
t r^ <y-°- 8 i n the shear thinning regime, in agreement 
with a recent experiment |36| . while r ~ 'j~ 2 - 5 in the 
shear thickening regime. The distinct shear rate scaling 
of the relaxation time together with the breakdown of the 
scaling collapse of short time MSD's are thus dynamical 
signatures of the shear thinning— thickening transition. 

In summary, in dense athermal shear flows, we observe 
robust signatures of the shear thinning— thickening tran- 
sition in the granular temperature, structure, and dy- 
namics. Preliminary results indicate that our findings 



may be general to various systems. The quantities sorted 
out through this work would be useful and accessible 
tools in experiments and systematic simulations to un- 
derstand shear thickening at high shear rates. 
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